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There is a renewed interest in dynamic characteristics of damaged aircraft both in order 
to assess survivability and to develop control laws to enhance survivability. This paper 
presents a set of flight dynamics equations of motion for a rigid body not necessarily 
referenced to the body’s center of mass. Such equations can be used when the body loses a 
portion of its mass and it is desired to track the motion of the body’s previous center of 
mass/reference frame now that the mass center has moved to a new position. Furthermore, 
results for equations presented in this paper and equations in standard aircraft simulations 
are compared for a scenario involving a generic transport aircraft configuration subject to 
wing damage. 


I. Introduction 

There is a renewed interest in control of aircraft that have sustained damage. This interest is driven by diverse 
factors such as advances in adaptive flight control, threats from shoulder-fired missiles, mid-air collisions, aging 
aircraft, and undetermined lifecycles of new composite transport aircraft. Illustrative examples are shown in figure 
1. In order to analyze the dynamics of damaged aircraft the dynamic equations of motion must properly reflect the 
underlying physics. This is not a novel problem from the dynamic perspective. Since the first payloads have been 
delivered from aircraft in flight, the equations of motion describing these dynamics have been derived from first 
principles. 1 With the renewed interest in the problem, the authors hope to provide a convenient reference for the 
dynamic equations of motion that properly account for non-collocation of the center of mass and the body axes 
reference point. 



(a) (b) 

Figure 1 (a) Air Transat A310 - 2005 - Rudder Loss; (b) Embraer Legacy 600 - Brazil - Wing Tip and 
Stabilator Loss, September 2006. 

To provide a clear exposition, the paper is organized as follows. First, rigid body equations of motion are 
developed that are not necessarily referenced to the body’s center of mass. These equations govern the motion not of 
the center of mass but some other selected fixed point on the body. One could, for example, place the equation of 
motion reference at the aerodynamic reference center and automatically obtain the velocity components required to 
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define angle of attack and sideslip angle, independent variables for the aerodynamic forces and moments. The 
proposed equations of motion are then applied to the problem of modeling large instantaneous shifts in center of 
mass of a rotating body that suddenly loses a portion of its mass. Benefits in simulation implementation are 
discussed. As an example, open loop dynamic responses of a generic transport aircraft with wing damage are 
compared for typical equations of motion used in the aircraft simulation and the equations presented here. 


II. Rigid Body Equations of Motion Referenced to an Arbitrary Fixed Point on the Body 

There are several approaches that can be used to develop the general equations of motion. The one selected here 
starts with Newton’s laws applied to a collection of particles defining the rigid body (any number of dynamics or 
physics books can serve as references, e.g. reference 2). In this paper, the rigid body equations of motion over a flat 
non-rotating earth are developed that are not necessarily referenced to the body’s center of mass. Such equations 
will be used in the next section when the body loses a portion of its mass and it is desired to track the motion of the 
body’s previous center of mass/reference frame now that the mass center has moved to a new position 


0) 



Figure 2. Arbitrary body referenced in inertial and body-centrie frames. 

A. Translational Acceleration 

Consider the rigid body in figure 2. Let the X-Y-Z reference frame be fixed inertially at point O and the x-y-z 
reference frame be attached to the body at any arbitrary point A. The body is free to move and has an instantaneous 
angular rate (O with respect to reference X-Y-Z. Let m i denote a particle of mass of the body located in the fixed 
inertial reference frame as 

r ,= r i+P, (!) 

where r 4 defines the position of point A with respect to O and p, defines the position of m with respect to point A. 
Furthermore, let G denote the center of mass located by p from A and r from O. From Varignon’s Theorem 2 , p 
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and r satisfy mp = X w,P, and mr ='£ J mr i , respectively, where m is the total mass of the body, m = X n h ■ 
From Newton’s Second Law, the translational motion of the body is governed by 

P=I»,(f,)«=-(F)„ (2) 

where X F is the sum of all external forces on m and (r ) is the acceleration of the center of mass observed from 
the fixed inertial reference frame. 

To rewrite (2) in terms of components measured relative to the body-fixed x-y-z reference at A, take the 
derivative of (1) with respect to X-Y-Z following the rule 

+ a> x [ ] so (r, ) xrz = (r 4 ) XYZ + (p,- ) XYZ becomes 

(ii)xyz =y A + P,+® X P, ( 3 ) 

where x A is (r A ) expressed in x-y-z coordinates. Here, p, and to are expressed in x-y-z components and p, is 
the velocity of m i as observed in the x-y-z frame. To promote economy of notation, ( ) will be dropped when 
the derivative is with respect to the x-y-z frame. Taking the derivative of (3) again with respect to X-Y-Z yields 


fc)^ = (^ + P* + ® x P/)J +«x(v 4 +p.+(oxp.) (4) 

or 

= v 4 + toxv 4 + p,. + (bxp. + 2-(a>xp.) + <ox(©xp.) (5) 

Let IF be expressed in components measured relative to the body fixed axis, then 

XF = 'Lm t (\ a +cox v 4 +p ( + to x p. + 2 • (to x p.) + to x (to x p.)) (6) 

XF = m(\ A + to x v,) + X w,p, + 2 • (to x X mp j ) + ti) x X m l p l +tox (to xXw ; p ( ) (7) 

XF = m(y A + tox v^) + mp + 2• (tox mp) + tox mp + tox (tox mp). ( 8 ) 


Since the body is rigid p = p = 0 , 

XF = m(\ A +iox v^) + tbx/Hp + tox(tox/np). (9) 

To formulate aircraft equations of motion for a non-center of mass reference, define the angular velocity, 
position of the center of mass, and the velocity at A in the body fixed axis x-y-z as 

to = Pi + Qj + Rk 
p = Axi + Ayj + Azk 
y A =U A \ + V A \ + W A V 

respectively. Here, the velocity components are subscripted with A to distinguish them from the usual velocity 
components considered at the center of mass. Of course these components are the velocity of the center of mass 
when A is the center of mass, p = 0 . Separating the force of gravity 
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W = -wgsin ©i + »/gcos©sin<l>j + ff/gcos©cos<Dk ) 


from the other external forces, the force driven aircraft equations of motion referenced to an arbitrary body fixed 
point A are 


ZF r = m(U A + QW A - RV A - ( Q 2 + R 2 )Ax + (QP - R)Ay + (RP + Q )Az + g sin 0) (10) 

ZFj, = m(V A +RU A -PW A + (PQ + R)Ax 

~(P 2 + R 2 )Ay + (QR -P)Az- geos©sin<D) 

ZF z =m(W A +PV A -QU A + (PR-Q)Ax 

+ (QR + P)Ay - (P 2 +(Z)Az-gcos©cos<D). 

Here, O and 0 are the standard roll and pitch Euler angles from the standard aircraft yaw-pitch-roll Euler rotation 
sequence defining the bodies orientation to the fixed inertial reference frame. 3 The angular acceleration terms in 
(10-12) and additional angular rate terms are due to the center of mass offset from the body fixed frame. Solution of 
the coupled translational and angular accelerations will be discussed after the equations governing the angular rates 
are formulated for the offset body fixed axes. 

B. Angular Acceleration 

Consider again the rigid body described in figure 2 where the X-Y-Z reference frame is fixed inertially at point 
O and x-y-z reference is attached to body at any arbitraiy point A. The body is free to move and has an 
instantaneous angular velocity, to, observed from the X-Y-Z reference frame. Again, let /// denote a particle of 
mass of the body located from the fixed frame as 


r ,= r . 4 +P, (!) 

where r 4 defines the position of point A with respect to O and p ( defines the position of m i with respect to point A. 
In this section, the momentum principle, relating the rate change of angular momentum to the external moments, is 
formulated for the moment of linear momenta taken about A. 

The absolute angular moment about A is defined as 

H, =Z(p,xm,v ; ) (13) 

where v ; = (r j , , z is the velocity of m. observed from the X-Y-Z reference frame. The derivative of H ( observed 
from the X-Y-Z reference frame is 


(H ,) xw = Z((p, ) xw x m i \ i ) + Z(P, x m, (v,.) xyz ) (14) 

where (p,) nz and (v,-) are derivatives observed from the X- Y-Z reference. From figure 2, (p, ) xrz = (F" ) xr/ 
and ( v ) = (f ) . Since (r) xv =0 and the sum of external moments about A is defined as 

ZM, =Z(Pi x m i (r ) ) , 2 (14) can be rewritten as 


ZM x=(^L + (^)^ xZw ' V '- 


( 15 ) 
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The velocity of the center of mass with respect to frame X-Y-Z, v, satisfies = Noting v = v 4 +(p) 

where both \ A (=(r A ) ii7 ) and (p j are observed from X-Y-Z, the desired momentum principle equation is 

zm -=KL +v ^ xw (pL- (16) 

If A is the center of mass, v equals \ A and (15) assumes the usual form, In order to express (16) 

in components measured in the body fixed frame x-y-z at A, first express (13) as 

H , = I(p, xot.(v 4 +p,. +©xp.)) 

where p, is the velocity of m j observed in frame x-y-z. Vectors p,, p,., to , and v A are all expressed in components 
measured in the body fixed x-y-z frame. Since the body is rigid, p, = 0 so 

H , = Y m i p, xv , + Ip, x m i (to x p . ) or 

H 4 = »ipx v 4 +Ioj (17) 

where I is the mass moment of inertia matrix of the rigid body. I is defined relative to the body fixed x-y-z at A and 
hence, is invariant with time, 1=0. Iw is a vector in the x-y-z frame arising from the integration version of 
Ip, xm,.(wxp.), jpx(coxp )dm and will be given below. Since and p are now expressed in components 
measured relative to frame x-y-z, the needed absolute derivatives measured from X-Y-Z in (16) are obtained from 

(18) 

and 

(pL=(pL +a,x p- (19) 

Substituting (17) into (18) yields 

=^(wp x v^+ I (»)J +wx(mpxv,+I(») 

-»<pxv 4 +mpx\ A + Id) + ci)x(wpx v^) + (0xl(a 


(h,) ve =Id) + ioxIci) + Hipxv i +wx(fflpx’vJ. (20) 

As before, derivatives with no noted referenced axes are assumed to be with respect to the body fixed x-y-z 
frame, i.e. ( ) . With regard to (19), the body is rigid, thus, the center of mass does not move relative to A. So 

p = 0 and (p) nz = tax p . Substituting this result and (20) into (16) yields 

EMj = Id) + 0) X 1(0 + »;p X V 4 +»kox(p X \ A ) + m\ A X (tox p). (21) 
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Redefining to not include the moment at A due to the gravity force acting through the center of mass yields 


ZMj = Id) + ft)xI© + mp x v A +m(ax(pxy A ) + m\ A x (ioxp)-p x W , (22) 

the desired moment equations of motion referenced to the arbitrary body fixed points. 

For the moment driven aircraft equations of motion, note I® is shorthand for 

I® = (I XX P - I xy Q - IJt)i + (~I xy P + r yy Q -I r _R )j + (~I XZ P-1 r Q + I :: R )k (23) 

with no assumptions of symmetry ; about the x-z plane (i.e., I * 0,1 ^ 0) . So (22) expands to 

P, = ij - K/2 ~ 4* + r xy PR - i x: PQ +(4 - 4 )QR + (R 2 - Q 2 )i„ 

+m((PV A -QU A + W A -gcos©cos®)Aj + (*4 - RU A -V A + gcos©sin<l))Az) (24) 

IM 4 = -I xy p + 1 yy Q - 4 R + 4 PQ - IJ2R + (4 -I ZZ )PR+{P 1 -R 1 )I XZ 


+m((QU A -PV A -W A +g cos © cos <5) Ax + (QW A - R V A + U A + g sin ©)Az) (25) 

IM, =-4^-42+4*+ 42*-4™+04-4)*2+(2 2 -* 2 )4 

+m((7?4 -P4 +4 - geos © sin 0)Ax + (i?4 -04 -gsin©)Ag) (26) 

As with the force equations, the moment equations contain coupling of both angular and translational accelerations 
due to the center of mass offset. 

III. Modeling Large Instantaneous Center of Mass Shift 

In the previous section, equations of motion were developed referenced to a point fixed on a rigid body, not 
necessarily the center of mass, with the center of mass fixed relative to this point. The problem of interest, however, 
is not one where the center of mass is fixed but one where at some point in time the center of mass instantaneously 
shifts to another location on the body. Consider a rigid body that at some time a piece of the body breaks away and 
motion of the remaining piece is of interest. Without loss of generality, assume that reference point A, discussed 
above, is attached to the body of interest. Examples include a fighter aircraft releasing its stores, a free falling space 
capsule separating from its launch abort tower, and an aircraft losing part of its wing. Momentum exchange due to 
mid-air collisions resulting in wing loss is not modeled here but can be included in the form of transient forces and 
moments. 

In this section, we contrast a more traditional center of mass (CM)-centric approach with one utilizing the 
equations developed above. The simulation modified for this study was of a subscale generic transport model that 
ran with fixed-time step, AT , Euler integration. Here, assume the center of mass shift occurs at some multiple of 
the time step, e.g. t = kAT . 

C. A. CM-Centric Approach 

The problem considered should not be confused with describing the motion of a rigid body with incremental 
changes in center of mass, e.g. fuel is burned. The general practice for incremental center of mass change is at each 
time step transfer the moments due to external forces referenced to some body fixed point to the new center of mass. 
The mass is modified along with the inertia tensor written with respect to the new center of mass, generally using the 
parallel axis theorem on constant mass subsystems. If the center of mass change is small, at each time step, the 
velocity at the center of mass of the previous time step is close enough to the velocity at the new center of mass so 
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the integration can run with that as the initial condition for the next integration. When the center of mass change is 
large, however, the velocity corresponding to the previous center of mass is appropriate as the next step initial 
condition for integration only when angular rate is near zero. If the body is experiencing any angular rate, as is 
common in aircraft maneuvering, such a strategy leads to poorly modeled transients since the initial condition does 
not reflect the velocity at the new center of mass. 

In the event of a large instantaneous change in center of mass, a carefully choreographed reset (initial 
condition change) of the integrators corresponding to v must occur. This reset involves 1) a trigger which detects 

the center of mass shift, 2) a calculation of v and v of the new center of mass location based on the center of mass 
shift and on previous values of the state and state rates of the previous center of mass location, 3) a prediction of the 
current v of the new center of mass location utilizing the results of 2), and 4) a reset of the integrator’s initial 
condition using this v of the new center of mass triggered by 1). The equations governing the motion of the center 
of mass can be obtained above by setting the center of mass offset from A to zero. Note, (A and to (as well as the 
Euler angles and rates) correspond to the rigid body, not specific to any fixed point on the body. As a result, the 
integrators on to need not be reset. 


In the simulation, the state [ \(kAT ), at(kAT) ] and the previous value of its rate [ v((£ - 1)AF), 
to((/c — 1) A A) ] are fed back to an auxiliary subsystem for defining angle of attack, sideslip, and their rates required 
as independent variables in the definition of aerodynamic forces and moments. Strictly speaking, these signals are 
used to determine the motion of the body-fixed aerodynamic reference point (defined in the wind tunnel test) as an 
appropriate definition of angle of attack and sideslip. In Euler integration, the held rates are the one-sided state 
derivatives at t = kAT approaching from left side. In this paper, both state and the previous rate will be considered 
current values and denoted as \ k cml and v* ml . Here subscripts end) and cm 1 denote its previous and new center of 
mass locations, respectively. Let and F*, define the position of the center of mass for each location, and 
AF = r k m , - F^, 1 define the change in center of mass with respect to its previous location. The center of mass 
location is generally determined in the mass and inertia subsystem of the simulation. The trigger detects a rise in 
norm of AF at t = kAT . To construct the state at cm 1 at / = kAT , previous values of the state and its rate defined 
at end) are used as follows. First the previous state and rate at cm 1 are constructed 




+ to 


1 x Ar 


v*“! 


= V *' 1 


+ (b* 1 x Ar . 


(27) 

(28) 


An Euler prediction provides an estimate of the current state at eg 1 at t = kAT . 


vl, £v*';!+tl:!A t 


(29) 


This estimate is used as an initial condition to reset the integrators of the translational state rates at t = kAT . The 
non-CM approach using the equations of the previous section avoids this reset problem entirely because reference 
point A is stationary on the body. 


B. Non-CM Approach 

The non-CM approach does not imply that the location of the center of mass is not needed. The location of the 
center of mass with respect to A, p, is required in the equations of motion. Let the forces and moments due to 

gravity be subsumed back into the external force and moment, and XM ,, respectively. Equations (10-12) & 
(24-26) can be concisely expressed in matrix notation 




otI 3 

D 

-1 ‘ 

" IF " 


mSl x 



cb 


D 

I 


. XM,_ 




(0 


( 30 ) 
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where 



0 

-mAz 

mAy 


‘ 0 

-R 

Q~ 

D = 

mAz 

0 

-mAx 

, = 

R 

0 

-p 


—mAy 

mAx 

0 


_-Q 

P 

0 

r 



“1 

r _ 



_ 1 


" 0 

-w A 

4 ' 


IXX 

-4 

- 4 ' 

W A 

0 

-U A 

,and I= 

-4 

4 

-4 

r v A 

U A 

0 


.-4 

-4 

4 


Here, the inertia matrix is expressed about fixed point reference A. If part of the wing of an aircraft breaks off at 
t = kAT with mass m tjp and center of mass located with respect to A at p tip , i.e. (A x lip ,Ay ljp ,Az lip 'j, and an inertia 
matrix defined about the tip’s center of mass, \ tip cg , the damaged aircraft (dac )’s mass, center of mass, and inertia 
matrix replace the nominal as 



m -> m dac = m 

~ m ,ip 


(32) 


P^Pdac= m P- m „pP,ip/ 

' m j„c 

(33) 


Iij ^ Iij,dac ^ij ^ij,tip 

i,j 

= x,y,z 

(34) 

^a,tip,> 

cm+ m tip( A jtp+ Ak ?ip) 1 


= x,y,z\i j &i ^ k 

(35) 

Iij,Up 

~ 1ij,tip,cm ^tip ^tip A/ tip 

i,j 

= x,y,z\i*j 

(36) 


Here, only the piece removed is subtracted from the inertia matrix. This is in contrast to the CM-centric approach 
which requires the inertias of both the undamaged vehicle and the missing tip to be transferred to the new center of 
mass and then subtracted. 


Perhaps the only disadvantage in using the non-CM equations of motion is the requirement that a 6x6 symmetric 
matrix be inverted to provide the state rates. Computationally this does not seem to be a burden for simulation 
environments such as SIMULINK. It should be pointed out, however, that the input to the subsystem defining angle 
of attack and sideslip, as well as other required variables are defined with respect to states and rates at the vehicle’s 
center of mass. This does not pose a problem since these values are easily obtained with 


WxAp, te 
4nl =^A+<i> kxA P d ac 


(37) 

(38) 


These calculations are performed upstream from the auxiliary and aerodynamic subsystems that need them. 
Furthermore the auxiliary system defining angle of attack and sideslip can be greatly simplified by choosing the 
aerodynamic reference center as A, since each is defined by the velocity of the aerodynamic reference. 

Other benefits of the non-CM approach include the following. Only one moment arm needs to be monitored— 
that corresponding to gravitational force applied at the center of mass. Assuming the wind tunnel data has been 
specified at some aerodynamic reference fixed to the vehicle, the moment arm from fixed point reference A to the 
aerodynamic reference is time invariant—ideally zero with the correct selection of A. Since the reference A is body 
fixed during the damage event, the velocity states at A provide the appropriate initial state for the next integration 
step. No reset of the integrator is needed! There is no need for the trigger logic. There is no need for predicting the 
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current state at the new center of mass from previous states and rates corresponding to the previous center of mass to 
obtain the appropriate initial condition reset in a timely manner. In the next section, an example will illustrate that 
the non-CM approach duplicates the results of the CM-centric approach. 


IV. Example: Dynamic Response of Generic Transport Model with Damage 

Wind tunnel test for several cases of damage was conducted in 2006 at NASA Langley 14x22 Subsonic 
Tunnel. Figure 3 illustrates different damage cases tested. Both static and dynamic aerodynamic data along with 
aerodynamic increments for damage are incorporated into the Generic Transport Model (GTM) simulation along 
with corresponding changes in aircraft weight and inertia. In this example, the aircraft is trimmed at level flight. At 
five seconds, the aircraft responds to a 1.5 deg aileron command to roll right for 5 seconds. At ten seconds, the roll 
rate is zeroed for the next 5 seconds, and then at 15 seconds, wings are leveled with a 1.5 deg. aileron command to 
roll left for 5 seconds. Prior to completing this maneuver at 19 seconds, 24% of the left wing’s semi-span breaks off 
with mass equal to 1.95% of the entire vehicle. 

The mass of the missing piece may not be representative of actually losing 24% of the left wing. The GTM 
is a 5.5% dynamically scaled version of a general full-scale transport. A CAD generated model of the missing piece 
incorporating ribs, stringers, spars, and skin for the full-scale aircraft was constructed. To obtain the subscale 
version, similitude 4 rules dictate that the mass and inertia of the full-scale vehicle be multiplied by (.055) 3 and 
(.055) 5 , respectively. To provide a 1.95% drop in the vehicle’s mass, a factor was applied to the density of the 
material used in the full scale model. This factor modified both the missing piece’s mass and inertias. For a vehicle 
with a 6.8488 ft wingspan, this center of mass shift was 0.041 inch forward to nose, 0.715 inch to the right wing, and 
0.052 inch up from the nominal position—not a lot of change, but enough to show both approaches yield the same 
result. To clearly illustrate differences in motion of the two reference points considered, reference A was arbitrarily 
selected to be 2.0 inches behind, 4.8 inches to the right, and 4.5 inches above the nominal center of mass. 

Figure 4 shows the results for the CM-centric approach in blue and non-CM approach in red for the 
translational velocities coming from the integrator on the left side and the velocities at the center of mass generated 
by both approaches on the right side. As expected the velocities at reference A differ from those at the center of 
mass, but the center-of-mass velocities calculated from the reference A states (i.e. equations (37-38)) overlay the 
CM-centric results. 


A closer look at the instant of failure, in figure 5, reveals that the Non-CM approach exactly captures the 
velocity of the CM at upset without the need to reset the integrator! The angular rates and Euler angles are shown in 
figure 6 also reveal duplicate results. 



Figure 3. Commercial transport damage research . 
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Figure 4. U, V, W from integrator for CM-centric (blue) and non-CM approach (red), left, and calculated at 
center of mass, right. 


lOof 13 


American Institute of Aeronautics and Astronautics 






















W body, fps V body, fps 


190.4 


» 190.2 

a 


o 

X 

3 190 


189.8- 1 - 

18.95 19 19.05 




Time, sec. 


190.4 


a 190.2 


O 

x 

6 

190 


189.8- 1 - 

18.95 19 19.05 




Time, sec. 


Figure 5. Closer view at failure of U, V, W from integrator for CM-centric (blue) and non-CM approach 
(red), left, and calculated at center of mass, right. 
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Figure 6. Angular rate and Euler angle comparison 


V. Conclusion 

The Non-CM approach for modeling a rigid body’s motion through an instantaneous shift in the body’s center 
of mass as a result of mass loss essentially duplicates the result of the more standard CM-centric approach to 
simulation. The Non-CM approach utilizes equations of motion expressed about a fixed point on the rigid body that 
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does not have to be the center of mass. If the body has any appreciable level of angular rate at the instant that a 
piece breaks off, the standard CM-centric approach must reset the integrator on v with the initial condition of v at 
the new center of mass location. This discontinuity is avoided entirely with the Non-CM approach, since reference 
A is fixed to the same point before and after the damage event. As a result the Non-CM approach does not require 
an integrator reset and all the numerical manipulation that involves: trigger detection, calculation of past motion 
states of the new center of mass prior to damage, and prediction of the current motion states at the new center of 
mass required to reset the integral’s initial condition. 

The Non-CM approach does require tracking of the body’s center of mass as does the CM-centric approach. 
The center of mass position appears in the Non-CM equations of motion and is required for solving the moment due 
to gravitational force. However, this is the only moment arm that changes. Both aerodynamic and propulsive forces 
and moments are generally written with respect to a fixed point on the vehicle. Since reference A is also fixed, the 
moment arms used to translate these forces and moments to reference A are time invariant and need only be 
calculated at the initiation of the simulation. The inertia matrix, defined about A, is easily modified without 
application of the parallel axis theorem to transfer the undamaged aircraft inertia to the new center of mass. Only the 
relative distance to the center of mass of the missing piece to A, its mass, and inertia of the missing piece defined 
about its own center of mass are needed to appropriately modify the inertia matrix to that of damaged aircraft about 
reference A. The new center of mass position is determined by the mass and center of mass locations of the 
undamaged aircraft and the missing wing tip. 

The Non-CM approach does come with a price. Instead of inverting a 3x3 inertia matrix, a 6x6 mass-inertia 
matrix of the coupled system must be inverted. Utilizing the matrix/vector formulation of Section III, SIMULINK, 
and computational environment available today, this does not appear to be a limiting factor. It took longer for this 
author to implement a correct integrator reset than to implement the equations of motion of the Non-CM approach. 
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